% fig3 分岔图(Bifurcation)
clear; clc; close all;

%% —— 参数设置 —— 
g11 = 1; g12 = 1; g21 = -2; g22 = 2;   
mu_vec = linspace(0,1,1000);          % μ 采样点
init_vals = [ 0.1,  0.1;              % 两组初始值
             -0.1, -0.1 ];
colors   = [1,0,1; 0.6,0.3,0];         % 对应品红、棕色

Niter   = 500;   % 总迭代步数
Ntrans  = 400;   % 丢弃热身步数

nMu = numel(mu_vec);
nIV = size(init_vals,1);

% —— 预分配存储 —— 
% y_collect{j,i} 存第 j 组初值在 mu_vec(i) 处的所有后热身 y 值（列向量）
y_collect = cell(nIV, nMu);

%% —— 并行计算 —— 
parfor i = 1:nMu
    mu = mu_vec(i);
    tmp = cell(nIV,1);
    for j = 1:nIV
        x = init_vals(j,1);
        y = init_vals(j,2);
        ys = zeros(Niter - Ntrans,1);
        idx = 0;
        for k = 1:Niter
            x_next = mu*x + g11*sin(x) + g12*sin(y);
            y_next = mu*y + g21*sin(x) + g22*sin(y);
            x = x_next;  y = y_next;
            if k > Ntrans
                idx = idx + 1;
                ys(idx) = y;
            end
        end
        tmp{j} = ys;
    end
    y_collect(:,i) = tmp; 
end

%% —— 绘图 —— 
figure; hold on;
for j = 1:nIV
    % 把所有 μ 点和对应的 y 值展平
    all_mu = repelem(mu_vec, Niter-Ntrans);
    all_y  = reshape(cell2mat(y_collect(j,:)), [], 1);
    plot(all_mu, all_y, '.', 'Color', colors(j,:), 'MarkerSize', 1);
end

%% —— 美化 —— 
xlabel('\mu', 'FontSize', 14);
ylabel('y',   'FontSize', 14);
title('fig3(a)：\mu–y 分岔图', 'FontSize', 16);
xlim([0,1]);
grid on;

% 图例
h = gobjects(nIV,1);
for j = 1:nIV
    h(j) = plot(NaN, NaN, '.', 'Color', colors(j,:), 'MarkerSize', 8);
end
legend(h, {'(0.1,0.1)', '(-0.1,-0.1)'}, 'Location','northeast');

hold off;
